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Closed hierarchies and non-equilibirum steady states of driven systems 
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We present a class of tractable non-equilibrium dynamical quantum systems which includes com¬ 
binations of injection, detection and extraction of particles interspersed by unitary evolution. We 
show how such operations generate a hierarchy of equations tying lower correlation functions with 
higher order ones. The hierarchy closes for particular choices of measurements and leads to a rich 
class of evolutions whose long time behavior can be simulated efficiently. In particular, we use the 
method to describe the dynamics of current generation through a generalized quantum exclusion 
process, and exhibit an explicit formula for the long time energy distribution in the limit of weak 
driving. 


Significant activity has been devoted to the study of 
quantum systems out of equilibrium, with a rapid in¬ 
crease in interest due to the relevance to experiments 
with ultra-cold atomic gases, whose coherent evolution 
may be effectively controlled and decoupled from dissi¬ 
pation to a heat bath m- Non equilibrium dynamics is 
typically studied in processes such as external driving, re¬ 
peated quantum measurements and quantum quenches. 
The fundamental question that arises in such cases is 
what is the long term behavior of the system: does it 
eventually reach a non-equilibirum steady state? What 
is the nature of such a state? 

In studying the aforementioned non-equilibrium situa¬ 
tions, some highly successful tools of equilibrium statis¬ 
tical physics, such as linear response theory, may easily 
fail. Thus, there is a need to develop new methods to 
deal with some of these problems. Here we focus on one 
such idea - that of establishing closed hierarchies in or¬ 
der to get tractable equations for correlation functions. 
Specifically, in many statistical mechanics problems, it 
is possible to make a systematic connection between the 
evolution of n body density functions with n -|- 1 density 
functions. A prime example for such a set of relations is 
the Bogoliubov-Born-Green-Kirkwood-Yvon (BBKGY) 
hierarchy, which is the essential structure leading to the 
Boltzmann equation. In the Boltzmann equation, single 
particle densities are tied to higher order correlation func¬ 
tions represented in the collision integral (see, e.g. [4]). 
In this letter, we describe the requirements on obtaining a 
hierarchy under general quantum operations on fermions. 
We then show how the hierarchy may be closed for a 
quantum system that is periodically evolved, detected, 
and injected with current. Finally, we use the idea to de¬ 
scribe dynamics of current buildup, and the energy dis¬ 
tribution in the long term non-equillibirum steady state. 

To begin the discussion, consider the most general evo¬ 
lution of a density matrix, describing unitary evolu¬ 
tion, measurements and interaction with the environ¬ 
ment. Written as 

p ^ C{p) = = 1 (1) 

This form ensures p remains a non-negative matrix, and 


the normalization condition on the Krauss operators Ai 
ensures that Trp = 1 is preserved under the evolution. 

In general, there is no simple relation between correlation 
functions computed in state p before and after the evolu¬ 
tion which necessitates working in an exponentially 
large Hilbert space and is therefore often un-tractable. 
Hierarchy structures have been used before in the context 
of Kossakowski-Lindblad evolution, which is a particular 
limit of 0 . For example, the steady state of a dissipative 
XX spin chain in the presence of driving and dissipation 
has been studied extensively [SHI]- Also, conditions for a 
closed hierarchy in the continuous time frame work where 
stated in [^. Here we concentrate on a discrete time 
framework, but also supply corresponding Kossakowski- 
Lindblad results in the end as a special limit. In other 
processes, the possibility of getting a closed equation for 
Kossakowski-Lindblad evolution of noise averaged expec¬ 
tation values was studied in [10], to explore the stability 
of fractional charges to noisy hopping processes. 

We utilize the power of this approach to study a non- 
equilibirum process of current generation, as schemati¬ 
cally depicted in Fig. (a). In this process, we connect 
site a to a lead, where a current is injected, and parti¬ 
cles are allowed to go out at site b (two choices for b are 
shown). The process is explicitly described by 

P — >U{{l-r)p + ra[alpaa + napna]+ (2) 

r(l - a)[abpal + (1 - nb)p{l - nb)])W, 


where checks for the presence of a 

fermion on the injection/extraction site, and U = 

evolution between attempts 
during a time interval r. Here r is the overall attempt 
rate, and a is the relative probability of injecting vs ex¬ 
tracting attempts. We show below that this process leads 
to a closed equation (13) for the two point function of the 
system, which can be then computed numerically. It is 
important to emphasize that the long time steady state 
reached by the system is not a thermal equilibrium state, 
in that the energy occupation is very different from a 
Fermi-Dirac distribution governed by h. 


For small r, we find a remarkable explicit formula for 
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the steady state distribution = {steady\al.ak\steady). 
Here k labels the eigenstates \k) of the single par¬ 
ticle hamiltonian hnm, h\k) = Ek\k). Let pa^k = 
\{a\k)\‘^,Ph^k = l(a|fc)P be overlaps of these states with 
the sites a, b. Then is a function of the ratio Pa,k/Pb,k' 


= 


A + B 


Pa,,k 

Pb,k 


(1 


a) -I- a 


Pa,k 

Pb,k 


( 3 ) 


The coefficients A, B are given below in Eq. ( [T^ . We em¬ 
phasize that this expression is valid for any system obey¬ 
ing the form ([^. The details of the distribution depend 
of sensitively on the choice of parameters. For illustra¬ 
tion, we consider hopping on a chain of length N, with 
the standard Hamiltonian Hhop = 

corresponding to Dirichlet boundary conditions. In this 
case Pa,k/Pb,k = sin^(^)/sin^(]^)- In Fig- [^we il¬ 
lustrate the result with N = 100, and injection ar a = 1. 
We evolve the system from an initial vacuum state at 
t = 0. The results for extraction at the final and penul¬ 
timate sites b = 100,99 respectively, show sensitivity to 
the choice of operation sites. The energy distribution is 
computed numerically at long times and is clearly seen 
to approach $ in the long time limit. We stress that 
once driving has stopped, the energy distribution $ will 
remain the stationary distribution under the subsequent 
free evolution. 

We now turn to establishing the framework for our pro¬ 
cesses. We consider a system of fermions on a lattice of 
N sites. In Q we take Krauss operators of the form 
where is an evolution under a non¬ 
interacting hamiltonian, and is a polynomial of order 
in fermion operators a\a. The evolution under C of 
a general correlation function. 


{a^-A^o-ik+i-a-rk+i) = Trpa|^...aJ^ai,_,i...ai,_,, (4) 

is given by 

{A^-Ak'^ik+i-^bk+i) —^ + ( 5 ) 

El, Tr pC/tmt [4^ , m^U^] 

where the normalization relation in 0 was used. 

The assumption that the t/j, are non interacting, means 
that U^aiUi, = Uu-,ijaj for some unitary matrix it,, G 
U{N). As a consequence the evolution of the k +1 corre¬ 
lation function Q, is related in ([^ to correlation func¬ 
tions of an order at most k + l + 2 establishing 

a hierarchy of equations. 

We emphasize that the resulting state may be arbitrar¬ 
ily complex. Indeed, even when starting with a non¬ 
interacting thermal state, p ^ exp{—hija\aj) and taking 
each a non interacting unitary, p evolves into a sum 
of exponentials of fermion bi-linears. Such a state can be 
used to approximate any interacting state whose deter¬ 
minant quantum Monte Carlo description does not suffer 
from a sign problem [In- 



Figure 1: (a) Fermions hopping on a chain. Particles are 
injected from the left and removed on the right with r = 
0.01,0 = 0.7, r = 0.1. (b) Energy occupation approach to 
Results for extraction site b — 100 (upper panel) and 
6 = 99 (lower panel). There are 300 iterations between suc¬ 
cessive curves. For reference a Fermi-Dirac distribution is also 
shown, (c) Evolution of the local density, (ajai), in real space 
(red/blue corresponds to high/low density). 
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Below, we list several fundamental operations under 
which the hierarchy closes at the two point function level, 
for Gij = {a\aj), inducing a map G —)■ IC{G). We start 
with the obvious one: 

(I) The non-interacting evolution £„(p) = UpW, as de¬ 
scribed above, induces a map 

Gij —>■ /Cu(G)ij = {v) Gu)ij ( 6 ) 

We augment the free evolution with the following types 
of operations acting on a single particle mode: particle 
detection, injection and extraction. Below, for simplicity 
of presentation we will associate the operation with the 
mode associated with site i. 

Denote Pi the matrix {Pi)ra = ^ir^ia the projection on 
site i, and P^ = 1 — Pi, we introduce: 

(II) Particle detection at site i: 


CD,i{p) = n^pni + {l-ni)p{l-ni) (7) 

where = ajai. The induced map on G is: 

ICdAG) = P^^GPi^ + P,GP, . (8) 

The process (II) may be viewed as a “decoherence” of the 
correlations G in between site i and the rest of the lattice. 
As a linear super-operator on matrices, the measurement 
ICD,i has a simple spectrum. It acts as identity on matri¬ 
ces which do not mix site i with the rest, hence the non¬ 
zero subspace of matrices has a dimension 1 -|- [dimP^)^. 
The complementary zero subspace is spanned by the off 
diagonal blocks, of dimensionality 2{dimP^). 

(Ill) Removal of a particle from site i is described by 

Gout,i{p) = + (1 - 'ni)p{l - Hi) (9) 

with the induced map on G: 

ICoutAG) = Pi^GPA ■ ( 10 ) 


present, the particle extraction map will generically drive 
G to 0, i.e. {ICulCout,iA 0 [21]. Similarly, adding par¬ 
ticles by injection {ICuICin,iA, with no extraction present, 
will result in Gij —>■ Sij, when n —>■ oo, which is the state 
where all sites are occupied. 

On the other hand the unitary evolution (J) and the de¬ 
tection process (//) preserve the average particle number, 
i.e. {J2i <^i^i} = Tr G remains constant under ICu,JCm- 
There are a myriad possible processes described by com¬ 
binations of the oprations {I — IV). Here we concentrate 
on current generation processes as described by Eq. ([^, 
involves operations I, III, IV resulting in the map: 

G^uH(l-r)G + r{aPi~GP^+ (13) 

(1 — a)PAGPA))u + raAPaU. 


This simple model allows for a substantial reduction of 
complexity from the full quantum problem of describing 
the evolution of p into an evolution equation for the two 
point function Gij, which can be tractable by either an¬ 
alytical or numerical methods. It is clear at this stage 
that we can access very interesting situations. 


To compute the eventual non-equilibrium steady state 
for (13) it is convenient to view the transformation on G 
from a point of view of a super-operator. Here the N x N 
matrix G is viewed as an N'^ dimensional vector, and the 
action of the evolution C on p translates in (13) into: 


G ^ AG -b 5, 


(14) 


where A is an x matrix, and g is the inhomoge¬ 
neous contribution due to the particle injection processes 
(12), and corresponding to the raAPaU term in (13). 

In general, whenever g = 0, the long time behavior will 
be determined as usual by the largest eigenvectors of A. 
However when g 0, the situation is somewhat differ¬ 
ent: Indeed, from Eq. (14), we see that when (1 — A) is 
invertible, there exists a unique stationarity G, that may 
be written in the form: 


As a super operator this simple map may be viewed as a 
projection on the space of matrices that do not have an 
(i,j) or {j,i) element for any j. 

(IV) Finally, this operation injects a particle at site i: 

Gin,i{p) = + "niPni ( 11 ) 

and induces the map 

IC,nAG)=A + PAGPA. ( 12 ) 

We note that in contrast with (/ — III), the injection 
ICin,i is an in-homogenuous transformation on matrices, 
a property which we use below to compute steady states. 

We can combine any of the site operations (H-IV) with 
the unitary evolutions (1) mixing the the addressed site 
i with the rest of the sites. When no particle injection is 


Gsteady = (1 - (15) 

If A — 1 is not invertible, i.e. there are steady states 
AGr = Gr, it means that the evolution u has an invariant 
subspace which does not include the sites a, b. In this 
case one has to work with a generalized inverse of (A — 
1). A steady solution can either not-exist, or be non¬ 
unique of the form Gsteady ^ Gr + {I — A)~^g. While 
inhomogenous equations are a common occurrence in the 
study of steady states in classical driven systems, they 
are used less in quantum processes, where evolution is 
unitary. A recent example of such a non-homogenous 
equation in a quantum context is the calculation of the 
expectation values of spin components in the steady state 
of a spin undergoing periodic laser pulses iniiisj. 

In the limit of r <C 1, we were able to solve exactly for 
the degenerate perturbation theory to lowest order in r. 
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obtaining for the energy distribution <i> the result (|^. 
While the result may be obtained by a systematic calcu¬ 
lation, it is rather lengthy, and instead here we simply 
state it as an anzats, which can be explicitly checked by 
the reader: taking G = J2k ^k\k){k\ +rD, where D is an 
off diagonal matrix, and plugging into (131 one can ver¬ 
ify that 4) is stationary under (131 to leading order in r. 
Some details of the derivation are supplied in the supple¬ 
mentary material. The coefficients are given by certain 
overlaps. Define: 


/i; = 2{apa,i + (1 - a)pb,i) ; Qst = , s,t e {a, 6} 

Then A, B in ([^ are given by 

A = ~ ■ B (16) 

where 


p — ^Qaa (1 Oi'jQbb T 0:(1 (Aji^QaaQbb Qab )) 

the form of the off-diagonal part D can also be obtained. 
We have verified the validity of the result numerically on 
numerous cases in addition to the one depicted in Fig. 
Sb). We see that to leading order, $ is independent 
of r. How can we understand this? Note that at r = 
0, there are infinitely many steady states (any G such 
that [G,h] = 0). However, when r fy 0, A stops being 
degenerate and it singles out a particular direction of 
breaking the degenerate space of matrices. 

The dependence of the dynamics on the initial condition 
is of interest by itself. While in Fig. we started the 
evolution from the vacuum state, in Fig. we describe 
such a process where the system is started off as the 
ground state of Hhop- The evolution happens in stages. 
In the initial stage of evolution we observe two shock 
wave fronts: one propagating with a region of reduced 
density from the right, collides with a front of enhanced 
density propagated from the left. It is interesting to note 
that the evolution is on a faster time scale than the speed 
of propagation of a wave-packet localized at a point by 
free evolution. In the context of classical non equilib¬ 
rium processes, shock waves have been described for the 
asymmetric exclusion process in e.g.|l4j (It is possible to 
use the present system also to describe such situations, 
however this will be done elsewhere). 

As the fronts collide the imbalance between the left and 
right sides of the chain starts to decrease. Finally, soliton 
like density packets of different velocities, are observed 
at longer time scales, and may be related to the soli¬ 
ton described in [15] in the context of the orthogonality 
catastrophe. It is interesting to note the injected par¬ 
ticles traveling from the left travel with faster velocities 
compared to their partners from the other side. 

In Fig. we show the average particle density n = 
N~^TrG. One of the interesting features observed is 
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Figure 2: Density depletion in a system where particles are 
extracted from the right at higher rate than injected on the 
left, here r = 1 and a = 0.3, initial state is the half filled 
ground state of Hhop- (a) Real space evolution of local density, 
(b) Evolution of space averaged density. 


a qualitative change in the slope of n(t) around 350 iter¬ 
ations. This change seems to correspond to the annihi¬ 
lation of the high density front coming from the left. To 
check this behavior, we consider, in Fig. [^the evolution 
when the initial stage is asymmetric itself: Here in the 
initial stage all sites i on the left, i < 100, are empty, 
while all sites on the right i > 100 are occupied. This 
state evolves through four fronts that collide and even¬ 
tually annihilate. Note that for coherent evolution from 
such an initial state, it has been shown that the front 
propagation has a scaling 1/fy [T^. In the context of 
evolution of magnetization in a spin chain the evolution 
of initial domain wall was studied in m 

Comparing the density evolution in Fig. [^and Fig. we 
see that there is a transient behavior associated with the 
different nature of the initial states, and their stages of 
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Figure 3: (a,b) Density evolution under the same dynamics as 
Fig-d however with a domain wall as an initial state. Here 
the depletion happens in two steps, but eventually reaches the 
same asymptotic value (h ~ 0.38) as for the initial uniform 
case. 


mined by the ratio between the effective tunneling prob¬ 
ability into energy k from site a compared to the effective 
tunneling rate of the state k through site b. 

Summary: We presented a class of non equilibrium quan¬ 
tum processes that correspond to closed hierarchies of 
evolution equations, and can thus be studied numer¬ 
ically efficiently. We have used this idea to explore 
non-equilibrium generation of currents and approach to 
steady states. We remark that the resulting states may 
also be viewed as Floquet states, and we have thus sup¬ 
plied a particular way of engineering such states, that 
may be of interest in the context of topological Flo¬ 
quet states [l8ll20j . Moreover, the energy distribution 
should be studied further: one can hope to test the re¬ 
sulting highly excited current carrying steady states in 
a variety of settings from cold atoms to mesoscopic sys¬ 
tems and spin chains. We emphasize that our result does 
not rely on integrability that has been used in studies of 
dissipative spin chains, and is available for periodically 
driven fermion systems that do not correspond to spin 
chains, including higher dimensional systems. 
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evolution. In Fig. there is a noticeable change in de¬ 
pletion rate around 100 and 300 iterations, the first kink 
corresponds the initial high density region on the right 
hitting the left side: at that point injection of particles 
becomes harder for a while and \dtn\ decreases until the 
density goes down enough on the left. The second kink 
is observed when the high density region is reflected back 
to the right: extracting particles on the right is then eas¬ 
ier and \dtn\ grows. At long times the density seems to 
decay asymptotically as 1 /t towards the non-equilibrium 
steady state density. 

Finally, we consider the Kossakowski-Lindblad limit, 
whose treatment is considerably simpler. Starting with: 


=-. (18) 

'~)aPa,k “t” ^bPb,k 

This last expression has a simple interpretation: the 
probability of occupying a given energy level is deter- 
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Next, we write: 

(Pa±$Pa±)„„ = ($ - Pa<^ " + Pa^Pa)nn = (24) 

‘^Pa.nn^n P 'P^lPa,n\^lPaAn 


and we write (23) as: 

0 = aPa.nn “ $n(aPa.nn + (1 “ a)Pb,nn) (25) 
+ S/($/ — 4 )„)(aPa,nlPa,ln + (1 — Oi) Pb,n\Pb,\n) 


where we also used: 

^lPa,n\Pa,\n — Pa,in 


(26) 


Supplementary Material 


In this section we derive the formulas for the en¬ 

ergy distribution $. 

We study the steady state equation associated with the 


process (13): 


G steady — (1 ^) G steady'^ P 

U^ra(Pa P Pa±G steady Pa P 

'a^r(l Cr) (P&yG^siearfy Pb_L )tt 


(19) 


At this point it is possible to argue that: 

— ^n){oiPa,n\Pa,ln + (1 — Oi) Pb.nlPbAn)\ (27) 
is small, giving us a first guess for the answer: 




(28) 


However, it is possible to do better and solve equation 


(23) exactly. 


To do so we write: 


where u = e . 

Below we label the eigenstates of ho by n, ho\n) = En\n), 
and would like to find the probability to find a state with 
energy En occupied in the steady state. This probability 
is given by = Tr{palan) = {n\G\n). 

For r = 0, all states where [G, h] = 0, are immediately 
invariant under time evolution. Therefore, in the limit of 
r <C 1 we look for an ansatz for the steady state Gsteady 
which is approximately diagonal. Let us write, in the 
energy basis, the ansatz: 


Gsteady = diag({$i,...})-brZ), (20) 


where = {n\Gsteady In) are the steady states occupa¬ 
tions, and D is an off-diagonal matrix in energy space. 
Eq. (19) becomes: 


$ + rD = (1 — r)<i) -b (1 — r)ru' Du P rau^ PaU +(21) 
rau"^{Pa±^Paj_)u P u'^r{l - a){Pbj_^Pb±)u P O(r^) 


We note that the zeroth order is eliminated and we wind 
up with: 

D = -^ + u^DuP (22) 

au'^PaU P au'^{Pa±^Pa±)u + u'l'(l - a)iPb±^Pb±)u 

Furthermore, note that both D,u^Du are off-diagonal 
in energy. Therefore we have a closed equation for the 
diagonal elements: 


Pa,ni = {n,a){aj) = f{n)f*{l) ; (29) 

Ph.ni = {n,b){bj) = g{n)g*{l) 


Going back to (23) we write it as: 


0 — Cy.Pa,Yin — 2$„(aPa,nn + (1 ~ Oi)Pb,nn) + (30) 

Ei^l{aPa,n\Pa.\n + (1 — Oi) Pb,n\Pb,\n) 


This equation becomes: 


0 = a|/np - 2^nia\fn\'^ + (1 - a)lg„l^) + 

^iMa\fn\WP{l-a)\gM^) 


(31) 


We rewrite the equation as a in-homogenous linear equa¬ 
tion: 

= aZpF + vt. (32) 

Here P is a unit vector defined by: 

= (33) 

Q is a diagonal matrix 

Qnm ~ ^nva's/JPn ? hn ~ 2(o; \fn\ + (7 )^34) 

and V can be written in the form 


0 — —+ aPa,nn + (23) 

Q:(Pa_L 4 )Pa_L)nn + (1 — 0 '){Pb±^Pb±)nn- 


Vnm = + (1 - a)l5np|5mp = (35) 

aZl\F){F\ + il-a)Z^a\G){G\. 
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The solution is given formally by: 


and we also need 


$ = 


(g2 _ ^ 

70.ZpF = aZpQr 


(36) 


Q2_y^^t^ 1-Q-iyQ 

Next, we define the unit vector \Fq) as 
\Fq) = Z-^Q-^\F). 

We see that: 

with the normalization 

^FQ = S„^; ||Fq|P = 1. 
Similarly we have 


-Q-^F. 


(37) 


(38) 


(39) 


= ^gqIGq); 


Using these, (36) is expressed as: 




$ = 

Zf Z^q 


)\Fq) 


-aZlZl^\FQ){FQ\-{l-a)ZlZl^\GQ){GQ\ 

In the next step we use the following relation: 

l+a|i?) (t!|+&|ii) (u| 1"^) 

l+a+b+aba-|(.,^)P) {(l + ^)l^) - b{u.v)\u)), 
which holds for normalized vectors ||m|| = I Ini I = 1. We 


(41) 


are not aware if the expression (41) appears in the litera¬ 


ture, but it can be verified explicitly by multiplying both 
sides by (1 + a|n)(n| + 6|u)(u|). 


We will use ( |4l| ) on (40), with |Fq), |Gq) playing the role 
of |m), |n). Thus, we take in (41): 


a —^ — aZ^ZpQ ; 


b —>■ —(1 — a)ZQZQQ, (42) 




ZqZfZfqZqq ' 


|gr.PlUI^ 


noting 


ZpZpQ = ; ZnZna = 


, -^G^GQ 


we have 


c = 


l9»I^IUI^ 




,4 /4, 

Sir ' ^ 


(44) 


(45) 


Using these expressions with (41) and (40) we find: 


^ _ aZpZpQ 

— 


UUt ((l-aZ^Z^Q|FQ)(FQ|-(l-a)zaz^Q|GQ)(GQ|)I^Q))’^ 

aZpZpQ _^_ 1 _ 

%/Ur l-aZ-i.Z-^-pZ-^Z-^^+aZ-j,Z-^pZ-^Z-^^{l-\cn ^ 

{(1 - /3Z2z2q)(„|Fq) +/3Z2z2QC*(n|GQ)} = 

(l-0Zlzl^)\f„\‘^+PZpZpQZGZGQC*\gr,\'^ 


Pn 1 aZj^ZpQ /3Z^ZQQ+aZpZpQ/3Z^ZQQ(l \c\^) 


(40) Explicitly: 


= 


so that: 

_ _a _ (1 —(1 —Q)Qbb)|/,iP + (l —a)Qba|g,iP (Ap,\ 

p„ l-aQaa-(l-a)Qbb + a(l-a)(Qaa<3bb-Qba^) ^ ' 

with 

^ ^ l//l^ ^ ^ 15/1"^ ^ ^ /.-N 

Gaa = - ; Gbb = - ( Gba = 1 

F F F 


which is the energy distribution formulas ( 3|16 ) in the 
paper. 





































